Monitoring-based post-earthquake damage assessment

0. Introduction

Near-real time damage assessment of the built infrastructure after earthquakes is a key enabler of rapid recovery and thus, resilience. Permanent monitoring of structures, for instance with acceleration sensors that have become increasingly available over the last years, offers data-sets before, during, and after earthquakes. Changes in the structural state due to nonlinearity and damage can therefore be picked up and correlated with the health of a structure.
Damage identification is commonly split into four levels (Rytter, 1993): 1) detecting the presence of damage; 2) locating the damage within the structure; 3) quantifying the amount of damage; and 4) assessing the influence of damage on the residual structural capacity. In the context of earthquake-damaged structures, the absence of detected damage (level 1) allows to declare buildings safe for occupance; the location of damage (level 2) is essential in planning and designing repair works; the severity of damage (level 3) is closely related to the widely applied damage classification in regional risk and loss-estimation applications; and the residual seismic capacity (level 4) governs the vulnerability and risk of a damaged structure with respect to aftershocks or subsequent earthquakes, for instance in a sequence.
While damage detection and localization do not necessarily require a physics-based model, damage quantification and, above all, residual capacity require fusing monitoring data with physics-based models.
In this demonstrator, the first three levels of damage identification are shown on a case study: shake-table tests of a half-scale four-story building with unreinforced masonry and reinforced concrete shear walls (Beyer et al., 2015). The file is organised as follows.
Table of Contents
0. Introduction 1. Case study description 2. Reference building state Compute the frequency ranges of interest 2.1 Derivation of the damage-sensitive features in healthy conditions Reference transmissibility Transmissibility Assurance Criterion Stiffness proxies Spectral energy distribution Modal identification 3. Data-driven damage detection 3.1 Mild earthquake (Test 2) Transmissibility and coherence Based on stiffness proxy Spectral energy distribution Modal identification 3.2 Moderately strong earthquake (Test 5) Transmissibility Based on stiffness proxy Spectral energy distribution Modal identification 3.3 Strong earthquake Transmissibility Based on stiffness proxy Spectral energy distribution Modal identification 3.4 Summary of data-driven damage detection 4. Data-driven Damage Localization 5. Damage quantification 5.1 Physics-based nonlinear simulation model References

1. Case study description

In order to demonstrate the capabilities of vibration-based monitoring for SHM and seismic condition assessment, in particular, the data from a shake-table test conduction by Tondelly et al. (2014) is used. This data corresponds to a half-scale unreinforced-masonry reinforced-concrete building. The test structure has been shaken by an earthquake with increasing amplitudes as shown in the figure below. Before each earthquake shaking a white-noise excitation has been introduced by the shake table.
The last earthquake led to the collapse (instability) of the structure and therefore, is not used here as monitoring is not required to identify near-collapse conditions. In other words, we are interested in forecasting damage prior to it becoming critical, and defining suited precursors for such a task.
clear; close all force; s = settings; s.matlab.editor.AllowFigureAnimation.TemporaryValue = 1;
load('WhiteNoiseData_PP.mat')
load('ConventionalData')
f_plotFigure0(WNM, DATASET)

2. Reference building state

Many damage detection (or anomaly detection) approaches involve the characterization of a healthy reference state of the structure. In this case, the white-noise test preceding the first earthquake excitation is used as reference state.
The reference signal is further divided into 3 segments, as shown in the figure below:
Given the limited dataset of white noise, some overlapping between the 3 segments is accepted. In real-world applications the healthy referenceprior to a seismic event can be assessed from sufficient ambiently collected data and thus, a much more comprehensive dataset would be available.
WNref = 1;
Signal = WNM.(['Test',num2str(WNref)]);
fsamp = 1024;
for Channel = 1:20
Signal(:, Channel) = Signal(:, Channel) - mean(Signal(:, Channel));
Signal(:, Channel) = f_butterworth(Signal(:, Channel), 256, 4, 45, 'low');
end
Ref_Signal0 = Signal(1*1024:21*1024, :);
Ref_Signal1 = Signal(16*1024:76*1024, :);
Ref_Signal2 = Signal(71*1024:191*1024, :);
f_plotfigure1(Signal, Ref_Signal0, Ref_Signal1, Ref_Signal2)
TransmPairs = [3 18; 6 18; 10 18; 14 18; 3 6; 6 10; 10 14];
TransmPairsRel = [1 5; 2 5; 3 5; 4 5; 1 2; 2 3; 3 4];
TransmPairNames = {'0-4', '1-4', '2-4', '3-4', '0-1', '1-2', '2-3'};

Compute the frequency ranges of interest

The transmissibility between two sensors gives an inditation of the frequencies, for which the signal is amplified. It is therefore linked to the stiffness of the structure above the reference sensor. The transmissibility is derived as follows:
,
where denotes the amplitude of the spectral density of the frequency range of the signal and i refers to the input signal, while o designates the output signal.
Earthquake signals are not stationary, therefore, wavelet decomposition is used here to transform the dynamic data from the time into the frequency domain.
In this first step, the first healthy dataset is used to derive the frequency ranges that contain structural modes. The quality of the damage-sensitive features that will be subsequently derived depends on the prominence of the transmissibility peaks, therefore only transmissibility peaks exceeding 3 are considered. The value of 3 (unitless as the transmissibility is the ration of two spectral amplitudes and characterizes the amplification of input signals) is chosen empirically, as values below get very noisy during the earthquake signal. The frequency bandwidth to calculate transmissibility is delimitated by the first crossing of 1 (indicating the end of the range that is amplified) and is slightly extended towards lower frequency values, given damage will reduce the stiffness of the structural system and thus, the natural frequency.
The transmissibility for four subsystems with the output sensor at the roof level and the input sensor (i) placed on the ground; (ii) placed on the first floor; (iii) placed on the second floor; and (iv) placed on the fourth floor are shown below. The identified peaks and the ranges to compute transmissibility-related quantities are also shown in the figure with vertical lines.
VPO = 48; TBW = 90; % Properties for wavelet decomposition
fanaRange = [1 40]; % Frequency bandwith that is considered
fb = cwtfilterbank('SignalLength',length(Ref_Signal0),'SamplingFrequency',fsamp,...
'FrequencyLimits', fanaRange, 'VoicesPerOctave',48, 'TimeBandWidth',90);
% Ranges are established from a classical PSD analysis
Channels = unique(TransmPairs);
nDOF = numel(Channels);
for j = 1:nDOF
[WTr.(['DOF',num2str(Channels(j))]), fr_WTr, CoI] = cwt(Ref_Signal0(:, Channels(j)), 'FilterBank', fb);
WTrAmp.(['DOF',num2str(Channels(j))]) = abs(WTr.(['DOF',num2str(Channels(j))]));
end
for j = 1:nDOF
WTr_Spectrum(:, Channels(j)) = mean(WTrAmp.(['DOF',num2str(Channels(j))])(:, find(CoI>1, 1, 'first'):find(CoI>1, 1, 'last')), 2);
end
Crit = [3, 2]; %critia for minimum peak height and range delimitation
fr_WTr_fl = flipud(fr_WTr); % flip the order as the wavelet decomposition produces vectors with decreasing frequency values (this is moslty for plotting reasons)
for TPdummy = 1:length(TransmPairs)
transm = WTr_Spectrum(:, TransmPairs(TPdummy,2))./WTr_Spectrum(:, TransmPairs(TPdummy,1));
transm = flipud(transm); %flip also the order of the spectrum
[pkvalues, pkindexes] = findpeaks(transm,'MinPeakHeight',Crit(1),'MinPeakDistance',15);
TrRef.nModes(TPdummy) = length(pkindexes); % for transmissibility
FSRef.nModes(TPdummy) = length(pkindexes); % for floor spectra
LMRef.nModes(TPdummy) = length(pkindexes); % for values related to the peak
for ModeDummy = 1:TrRef.nModes(TPdummy)
Ind1 = find(transm(1:pkindexes(ModeDummy))<Crit(2),1,'last');
TrRef.(['Pair',num2str(TPdummy)]).fRange(ModeDummy, 1) = fr_WTr_fl(Ind1) - 2;
Ind2 = pkindexes(ModeDummy) + find(transm(pkindexes(ModeDummy):length(transm))<Crit(2),1,'first');
TrRef.(['Pair',num2str(TPdummy)]).fRange(ModeDummy, 2) = fr_WTr_fl(Ind2);
LMRef.(['Pair',num2str(TPdummy)]).Peak(ModeDummy) = fr_WTr_fl(pkindexes(ModeDummy));
FSRef.(['Pair',num2str(TPdummy)]).fRange(ModeDummy, 1) = 0.5*TrRef.(['Pair',num2str(TPdummy)]).fRange(ModeDummy, 1) + 0.5*fr_WTr_fl(pkindexes(ModeDummy));
FSRef.(['Pair',num2str(TPdummy)]).fRange(ModeDummy, 2) = 0.5*TrRef.(['Pair',num2str(TPdummy)]).fRange(ModeDummy, 2) + 0.5*fr_WTr_fl(pkindexes(ModeDummy));
end
end
f_plotFigure2(WTr_Spectrum, TransmPairs, TransmPairNames, fr_WTr_fl, TrRef, LMRef)
As can be seen in the figure above, only the transmissibility between the ground and the roof (0-4) presents two modes in the frequency range 1-40 Hz that exceed the threshold value of 3. The red lines delimitate the frequency range that will be used to assess the changes in transmissibility. The dotted line indicates the frequencies, for which the amplification is highest. For instance, the modes of are 7.4Hz and 25.2 Hz.

2.1 Derivation of the damage-sensitive features in healthy conditions

Reference transmissibility

Based on the ranges defined in the previous part, the reference transmissibility is computed with the second reference dataset. The transmissibility of the seven considered sensor pairs is shown below. A Savitzky-Golay filter is used to smooth the wavelet transform.
fb = cwtfilterbank('SignalLength',length(Ref_Signal1),'SamplingFrequency',fsamp,...
'FrequencyLimits', fanaRange, 'VoicesPerOctave',48, 'TimeBandWidth',90);
for j = 1: nDOF
[WTr.(['DOF',num2str(Channels(j))]), fr_WTr, CoI] = cwt(Ref_Signal1(:, Channels(j)), 'FilterBank', fb);
DoNotDiscard = find(~(CoI>1));
WTrAmp.(['DOF',num2str(Channels(j))]) = abs(WTr.(['DOF',num2str(Channels(j))])(:, DoNotDiscard));
end
fr_WTr = flipud(fr_WTr); % flip the wavelet based spectra to have them with increasing frequencies
for j = 1:length(TransmPairs)
trans = mean(WTrAmp.(['DOF',num2str(TransmPairs(j,2))]),2) ./ mean(WTrAmp.(['DOF',num2str(TransmPairs(j,1))]), 2);
trans = flipud(trans);
TransmRef.(['Pair',num2str(j)]) = sgolayfilt(trans, 3, 9); % apply a filter for smoothing
end
f_plotFigure7(TrRef, TransmPairs, fr_WTr_fl, TransmPairNames, TransmRef)
As can be seen, the sensor pairs with adjacent floors (0-1, 1-2, 2-3) have a less pronounced transmissibility than the other sensor pairs.

Transmissibility Assurance Criterion

A first damage-sensitive feature (DSF) that can be used to detect the onset of damage is the Transmissibility Assurance Criterion (TAC) that returns the colinearity of the reference transmissibility (computed in the previous section) and the transmissibility of a new - potentially damaged - data window (Zhou et al, 2015). The formula of the TAC is shown below:
,
where is the transmissibility between sensor i and o; evaluated for the frequency range , characterizing the mode. The superscripts r and e denote the reference state and the new state to be evaluated, respectively.
While the wavelet decomposition theoretically provides a frequency-domain spectrum for each measurement sample (within the cone of influence); the frequency-domain spectrum is averaged over time windows (a duration of 5 seconds is chosen here) to reduce the influence of sensor noise and ground-motion variability.
In the following, the TAC is computed for the third healthy reference dataset and the corresponding distribution of the first four input-output pairs is shown below. The TAC equals 0 when the 2 transmissibility spectra are perfectly colinear and 1 when they are perpendicular.
% Chose parameters
WL = 5*fsamp; % Window length
OL = 2*fsamp; % Overlap between 2 windows
% % % % % % % % % % % %
clearvars WTr WTrAmp fb
fb = cwtfilterbank('SignalLength',length(Ref_Signal2),'SamplingFrequency',fsamp,...
'FrequencyLimits', fanaRange, 'VoicesPerOctave',48, 'TimeBandWidth',90);
for j = 1: nDOF
[WTr.(['DOF',num2str(Channels(j))]), fr_WTr, CoI] = cwt(Ref_Signal2(:, Channels(j)), 'FilterBank', fb);
WTrAmp.(['DOF',num2str(Channels(j))]) = abs(WTr.(['DOF',num2str(Channels(j))]));
end
for j = 1: nDOF
WTr.(['DOF',num2str(Channels(j))])(:, CoI>1) = [];
WTrAmp.(['DOF',num2str(Channels(j))])(:, CoI>1) = [];
end
fr_WTr_fl = flipud(fr_WTr);
%% Reference values and distributions for TAC
nWinds = floor(length(WTrAmp.DOF3)/WL)-1;
REF_FRAC.Mode1 = zeros(nWinds, length(TransmPairs) );
REF_FRAC.Mode2 = zeros(nWinds, length(TransmPairs) );
StartP = 1; EndP = WL; WindowDummy = 1;
while (EndP < length(WTr.(['DOF',num2str(Channels(j))])) )
% Reference value of TAC
for Pairdummy = 1:length(TransmPairs)
trans = sgolayfilt(WTrAmp.(['DOF',num2str(TransmPairs(Pairdummy,2))])(:, StartP: EndP), 3, 9) ...
./ sgolayfilt(WTrAmp.(['DOF',num2str(TransmPairs(Pairdummy,1))])(:, StartP: EndP), 3, 9);
trans = flipud(trans);
[MScoh, fmsc] = mscohere(Ref_Signal2(StartP: EndP, TransmPairs(Pairdummy, 1)), Ref_Signal2(StartP: EndP, TransmPairs(Pairdummy, 2)), 2*fsamp, 1*fsamp, 10*fsamp, fsamp);
for ModeDummy = 1:TrRef.nModes(Pairdummy)
Index1 = f_findmatch(fr_WTr_fl, TrRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1));
Index2 = f_findmatch(fr_WTr_fl, TrRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2));
REF_FRAC.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy) = 1- f_MAC(...
TransmRef.(['Pair',num2str(Pairdummy)])(Index1:Index2), ...
mean(trans(Index1:Index2, :), 2));
Index1 = f_findmatch(fmsc, TrRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1)); Index2 = f_findmatch(fmsc, TrRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2));
REF_MSC.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy) = 1 - mean(MScoh(Index1:Index2));
end
end
StartP = StartP + (WL -OL); EndP = EndP + (WL- OL); WindowDummy = WindowDummy + 1;
end
f_plotFigure5(TransmPairs, TrRef, REF_FRAC, TransmPairNames)
It can be seen that, despite the white noise excitation and absence of damage some variability in the TAC can be observed. This variability results from changes in the excitation between the (short) time windows and sensor noise.

Stiffness proxies

A second damage-sensitive feature is constructed as an approximator of the stiffness of the structural system. The proxy is computed as follows:
,
where is the acceleration measured at instance t at location dof. The stiff slab and the uni-axial shaking that characterize the test specimen result in the choice of one measurement point per floor and thus, dof coincides with the floor number, from 0 at the base of the building to 4 at the level of the roof. is the estimated displacement, which is obtained through numerical integration.
In a similar manner to the TAC, the value of the stiffness proxy (KPRX) is averaged over a time-window. In order to reduce the influence of signal noise on both the numerical integration that is performed to approximate displacement and the quantity that is correlated to force (summation of accelerations), the signal is reconstructed in a limited bandwidth around the structural modes. Therefore, the inverse wavelet transformation is used.
The KPRX for the first three sensor pairs is shown below. As can be seen, the proxy is rather stable when considering the entire system (), but the quality decreases rapidly when the size of the subsystem is reduced. Indeed, the linear fit - computed as the root-mean-square error of the linear prediction with the approximated stiffness (normalized with respect to the amplitude of acceleration) - is acceptable for the sensor pair (0-4) but increases rapidly, as can be seen on the bottom right plot.
Time = (0:1/fsamp:(length(WTr.DOF3)-1)/fsamp);
for Pairdummy = 1:length(TransmPairs)
for ModeDummy = 1:LMRef.nModes(Pairdummy)
sigtop = icwt(sgolayfilt(WTr.(['DOF',num2str(TransmPairs(Pairdummy,2))]), 3, 9), ...
fr_WTr, [FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1), FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2)]);
sigbot = icwt(sgolayfilt(WTr.(['DOF',num2str(TransmPairs(Pairdummy,1))]), 3, 9), ...
fr_WTr, [FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1), FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2)]);
DummyCounter = 1; IndexPairLow = find(Channels == TransmPairs(Pairdummy, 1));
for dummyLoop = IndexPairLow : nDOF
sigacc2(:, DummyCounter) = icwt(sgolayfilt(WTr.(['DOF',num2str(Channels(dummyLoop))]), 3, 9), ...
fr_WTr, [FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1), FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2)]);
DummyCounter = DummyCounter + 1;
end
sigaccFI = sum(sigacc2, 2); %Proxy of force (normalized by mass - which is considered unknown)
TopDisp = inte_frequency([Time', sigtop'], FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1), FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2), 2);
BaseDisp = inte_frequency([Time', sigbot'], FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1), FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2), 2);
RelDis = TopDisp - BaseDisp; %Proxy of relative displacement
[~, FitQuality.(['Mode',num2str(ModeDummy)])(:,Pairdummy) , Kproxy.(['Mode',num2str(ModeDummy)])(:,Pairdummy) , TimeVals, KprxOverall.(['Mode',num2str(ModeDummy)])(Pairdummy) ] = f_DSF_LC_batch(-1.*sigaccFI, RelDis', Time', fsamp, [1 35], WL, WL-OL);
end
end
f_plotFigure8(Kproxy, FitQuality)

Spectral energy distribution

The distribution of energy in the frequency domain can be used to compute another DSF. Based on the DSF proposed by Noh et al (2011), the ratio of the frequency of a fundamental mode of vibration is compared with the spread of energy in the surrounding frequency bandwidth:
,
where ψ is the energy at discrete wavelet coefficients. m denotes the wavelet coefficient of the m-th vibration mode and is the energy at n discrete wavelet coefficients that are used as reference spread around the vibration mode. These coefficients are chosen equidistanctly in the range [0.66 , ], where is the (pseudo-)frequency of the m-th mode. To avoid an excessive correlation with the frequency content of the ground motion, this DSF is reformulated in relative terms as the ratio between the output and the input channels.The probability density function (PDF) of the first three sensor pairs is shown in the figure below.
for j = 1:length(TransmPairs)
for ModeDummy = 1:FSRef.nModes(j)
[~, M1EthOUT, ~] = f_WTindex_M1E_batch(WTr.(['DOF',num2str(TransmPairs(j,2))]), fr_WTr,LMRef.(['Pair',num2str(j)]).Peak(ModeDummy), WL, Time, (WL-OL));
[~, M1EthIN, tvals] = f_WTindex_M1E_batch(WTr.(['DOF',num2str(TransmPairs(j,1))]), fr_WTr,LMRef.(['Pair',num2str(j)]).Peak(ModeDummy), WL, Time, (WL-OL));
REF_MER.Mean(j, ModeDummy) = mean(M1EthOUT ./ M1EthIN);
REF_MER.StDev(j, ModeDummy) = std(M1EthOUT ./ M1EthIN);
REF_MER.TimeHist.(['Mode',num2str(ModeDummy)])(:, j) = M1EthOUT ./ M1EthIN;
end
end
f_plotFigure10(REF_MER)

Modal identification

Modal properties, such as natural frequencies, consist another set of potential DSFs. Given the non-stationary nature of earthquake excitation, modal properties are derived from ambient vibrations not only for the reference state, but also later on for the potentially damaged state. In addition, the window length is extended for the modal properties to improve the stability in the identified features, which are the first two natural frequencies in the direction of shaking and the mode-shape curvatures of the corresponding modes.
The frequency-domain decomposition method is used here, as the controlled input facilitates automated peak-picking. Any other output-only modal analysis technique can be used to identify the modal properties.
The two identified frequencies and the corresponding mode shapes are shown in the figure below.
% Natural frequency and mode-shape curvature
[CPSD, fff] = f_CPSD_wind(Ref_Signal1(:, Channels), fsamp, 10*fsamp, 4*fsamp, 'hamming', 10*fsamp, 'c');
[Umat, Svec] = f_svd4fdd(CPSD, 1, 'act');
Svec = sgolayfilt(Svec, 3, 9);
EndIndex = f_findmatch(fff, 40);
[~, Peakindex] = findpeaks(Svec(1:EndIndex), "NPeaks", 2, "SortStr", 'descend', "MinPeakDistance", 20);
fnat_id = fff(sort(Peakindex, 'ascend'));
MScompl = Umat(sort(Peakindex, 'ascend'), :, 1)';
% get the relative distance from the ground point
MSrel = zeros(5,2);
for j = 2:5
MSrel(j,1) = norm(MScompl(1,1)-MScompl(j,1)); MSrel(j,2) = norm(MScompl(1,2)-MScompl(j,2));
end
% get the relative angle from the ground point
MSangle = angle(MScompl);
MSangle = MSangle - repmat(MSangle(1,:), 5, 1);
% construct the mode shape
Msha_id = MSrel .* sign(cos(MSangle));
Msha_id = Msha_id ./ repmat(max(abs(Msha_id)), 5, 1);
f_plotFigure11(fff, EndIndex, Svec, Msha_id, Peakindex)
Based on the identified modes (, peak picking in reduced frequency ranges is performed for shorter time windows to estimate reference distributions for both: the natural frequency and the mode-shape curvature. The mode-shape curvature is approximated by numerical derivation using Newton's formula:
, where x(i) is the modal displacement at the i-th degree of freedom and h is the height between two monitored degrees of freedom.The PDF of the two DSFs is shown in the figure below.
% get the mode-shape curvature
WL = 25*fsamp; OL = 20*fsamp;
StartP = 1; EndP = WL; WindowDummy = 1;
while EndP <= length(Ref_Signal2)
[CPSD, fff] = f_CPSD_wind(Ref_Signal2(StartP:EndP, Channels), fsamp, 10*fsamp, 7.5*fsamp, 'hamming', 10*fsamp, 'c');
[Umat, Svec] = f_svd4fdd(CPSD, 1, 'act');
Svec = sgolayfilt(Svec, 3, 9);
EndIndex = f_findmatch(fff, 40);
% as we have predefined natural frequencies, we conduct a guided search in limited bands around these frequencies
Peakindexes = zeros(1,2);
for ModeDummy = 1: 2
ind1 = f_findmatch(fff, fnat_id(ModeDummy)*0.66);
ind2 = f_findmatch(fff, fnat_id(ModeDummy)*1.25);
[~, PeakIndexes(ModeDummy)] = max(Svec(ind1: ind2));
PeakIndexes(ModeDummy) = PeakIndexes(ModeDummy) + ind1 -1;
end
Ref_fnat.Ref0.fnat_id(WindowDummy, :) = fff(PeakIndexes);
MScompl = Umat(sort(Peakindex, 'ascend'), :, 1)';
% get the relative distance from the ground point
for j = 2:5
MSrel(j,1) = norm(MScompl(1,1)-MScompl(j,1)); MSrel(j,2) = norm(MScompl(1,2)-MScompl(j,2));
end
% get the relative angle from the ground point
MSangle = angle(MScompl);
MSangle = MSangle - repmat(MSangle(1,:), 5, 1);
% construct the mode shape
Msha_sample = MSrel .* sign(cos(MSangle));
Msha_sample = Msha_sample ./ repmat(Msha_sample(5,:), 5, 1);
for ModeDummy = 1:2
for ChanDummy = 2:numel(unique(TransmPairs))-1
MSCref.(['Mode',num2str(ModeDummy)])(WindowDummy, ChanDummy-1) = (Msha_sample(ChanDummy-1, ModeDummy) - 2*Msha_sample(ChanDummy, ModeDummy) + Msha_sample(ChanDummy+1, ModeDummy)) / (1.475*1.475);
end
MSCref.(['Mode',num2str(ModeDummy)])(WindowDummy, numel(unique(TransmPairs))-1) = (Msha_sample(ChanDummy-1, ModeDummy) - Msha_sample(ChanDummy, ModeDummy) ) / (1.475*1.475);
end
StartP = StartP + (WL - OL); EndP = EndP + (WL - OL); WindowDummy = WindowDummy + 1;
end
f_plotFigure12(Ref_fnat, MSCref)
When using time window, some variability in the idenfied frequencies can be observed. In addition, the mode-shape curvature - despite the white noise conditions - show some variability.
clearvars sigaccFI sigacc2 WTr WTrAmp; save('ReferenceSignal0', 'MSCref', 'TransmRef', 'REF_FRAC', 'REF_MER', 'REF_MSC', 'Ref_fnat','FitQuality', 'Kproxy')

3. Data-driven damage detection

The next step consists in calculating the damage-sensitive features in near-real time when the earthquake occurs and the structural response is measured in the building.
To avoid overloading this demonstrator, three earthquakes test runs are used:
  1. Test #2 - leading to very minor damage (DG1 on the EMS scale)
  2. Test #5 - leading to moderate damage (DG2 in the EMS scale)
  3. Test #8 - leading to extensive damage (DG3 in the EMS scale)
The full analysis with additional details can be found in Reuland et al. (2022).

3.1 Mild earthquake (Test 2)

In real conditions, the signals characterizing the earthquake would be part of a permanent datastream. However, some DSFs and intensity measures (IMs) that may be of interest cannot be calculated for short moving time windows, but require a data window that covers the entire earthquake (without too much white noise prior to the earthquake and after the earthquake. The time values delimitating such a window are estimated using a STA-LTA algorithm that compares a short-term average (STA) and a long-term average (LTA) and allows to detect strong motions and delimitate them in time.
In the figure below, monitored accelerations during the second test are shown and the window that will be used for analysis are highlighted.
SigData = DATASET.Test2;
TimeVect = (1:length(SigData))./fsamp;
for ChanDummy = 1:nDOF
SigData(:, Channels(ChanDummy)) = f_butterworth(SigData(:, Channels(ChanDummy)), fsamp, 4, 35, 'low');
end
STA = 0.25; LTA = 5;
T1 = movavg(TimeVect', 'simple', floor(LTA*fsamp), 'fill');
T2 = movavg(TimeVect', 'simple', floor(STA*fsamp), 'fill');
MA1 = rms(SigData(:,1)) + movavg(SigData(:,1), 'simple', floor(LTA*fsamp), 'fill');
MA2 = rms(SigData(:,1)) + movavg(SigData(:,1), 'simple', floor(STA*fsamp), 'fill');
Ratio = abs(MA2) ./ abs(MA1);
Ratio = movmean(Ratio, 5);
Ind1 = min(find(Ratio < 0.7, 1, 'first'), find(Ratio > 1.5, 1, 'first')) - fsamp;
Ind2 = max(find(Ratio < 0.7, 1, 'last'), find(Ratio > 1.5, 1, 'last')) + 5*fsamp;
Ind2 = max(Ind2, Ind1+8*fsamp); % Impose a minimum length of 8 seconds
AnaRange = (Ind1 : Ind2);
f_plotFigure6(SigData,Channels,TimeVect,Ind1,Ind2,fsamp)
Once the signal delimitated, a wavelet decomposition is performed to analyse the signal in the time-frequency domain.
Sig2Analyse = SigData(Ind1: Ind2, :);
fb = cwtfilterbank('SignalLength',length(SigData),'SamplingFrequency',fsamp,...
'FrequencyLimits', fanaRange, 'VoicesPerOctave',48, 'TimeBandWidth',90);
for j = 1: nDOF
[wtrInstance, fr_WTr, ~] = cwt(SigData(:, Channels(j)), 'FilterBank', fb);
WTr.(['DOF',num2str(Channels(j))]) = wtrInstance(:, Ind1:Ind2);
WTrAmp.(['DOF',num2str(Channels(j))]) = abs(wtrInstance(:, Ind1:Ind2));
end
fr_WTr_flud = flipud(fr_WTr);

Transmissibility

The DSFs that have been described in the previous section are computed on the recorded resopnse during the seismic event. Unlike ambient signals, earthquakes are short signals, therefore the overlapping between time windows is increased to provide a sufficient amount of data points. The time increment between two analysis windows is chosen as 0.5 seconds, leading to two datapoints per second.
In addition to the DSF, the Mahalonobis distance to the healthy reference distribution is computed for each DSF. This distance allows to measure the discrepency between the new datapoint and the reference distribution (Mahalonobis, 1936).
WL = 5*fsamp; OL = 4.5*fsamp;
StartP = 1; EndP = WL; WindowDummy = 1;
while EndP < length(Sig2Analyse)
for Pairdummy = 1:length(TransmPairs)
trans = sgolayfilt(WTrAmp.(['DOF',num2str(TransmPairs(Pairdummy,2))])(:, StartP:EndP), 3, 9) ...
./ sgolayfilt(WTrAmp.(['DOF',num2str(TransmPairs(Pairdummy,1))])(:, StartP:EndP), 3, 9);
RefRatio = 1 ./ sgolayfilt(WTrAmp.(['DOF',num2str(TransmPairs(Pairdummy,1))])(:, StartP:EndP), 3, 9); % this is used to avoid strong influence of the input spectrum
trans = flipud(trans);
[MScoh, fmsc] = mscohere(Sig2Analyse(StartP:EndP, TransmPairs(Pairdummy, 1)), Sig2Analyse(StartP:EndP, TransmPairs(Pairdummy, 2)), 2*fsamp, 1*fsamp, 10*fsamp, fsamp);
for ModeDummy = 1:TrRef.nModes(Pairdummy)
Index1 = f_findmatch(fr_WTr_flud, TrRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1));
Index2 = f_findmatch(fr_WTr_flud, TrRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2));
DSF_TAC.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy) = 1- f_MAC(...
TransmRef.(['Pair',num2str(Pairdummy)])(Index1:Index2), ...
mean(trans(Index1:Index2, :), 2));
RefWeight.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy) = f_MAC(...
mean(trans(Index1:Index2, :), 2), ...
mean(RefRatio(Index1:Index2,:),2));
DSF_TAC.mahadi.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy) = mahal(DSF_TAC.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy), REF_FRAC.Mode1(:, Pairdummy));
Index1 = f_findmatch(fmsc, TrRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1)); Index2 = f_findmatch(fmsc, TrRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2));
DSF_MSC.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy) = 1 - mean(MScoh(Index1:Index2));
DSF_MSC.mahadi.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy) = mahal(DSF_MSC.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy), REF_MSC.(['Mode',num2str(ModeDummy)])(:,Pairdummy));
end
end
StartP = StartP + (WL-OL); EndP = EndP + (WL-OL); WindowDummy = WindowDummy + 1;
end
Mahal.TAC.EQK2 = DSF_TAC.mahadi;
In addition to the DSF samples calculated for short time windows of 5 seconds, single datapoints that average the transmissibility assurance criterion over the entire earthquake signal are computed. To limit the influence of the transient and non-white excitation, the TAC instances are weighted with a criterion based on the colinearity of the transmissibility and the ration 1/InputSpectrum. Finally, a multi-dimensional Mahalonobis distance is calculated for all 7 sensor combinations.
for Pairdummy = 1:length(TransmPairs)
DSF_TAC.OverallWeighted.(['Mode',num2str(ModeDummy)])(Pairdummy) = ...
sum(DSF_TAC.(['Mode',num2str(ModeDummy)])(:,Pairdummy) .* ...
(1-RefWeight.(['Mode',num2str(ModeDummy)])(:,Pairdummy))) / ...
sum((1-RefWeight.(['Mode',num2str(ModeDummy)])(:,Pairdummy)));
end
DSF_TAC_OverallMaha = mahal(DSF_TAC.OverallWeighted.Mode1, REF_FRAC.Mode1);
In the figure below, the datasamples computed for short time-windows during the earthquake are displayed and compared with the 95-percentile confidence interval of the healthy reference data (in green). As can be seen, the datapoints fall outside the reference range indicating:
The fact that the DSF falls into the healthy confidence range when the amplitude of shaking decrease indicates that some minor cracks may form during the shaking, but the close again when the excitation is reduced. Therefore, even though some minor damage is detected, no residual effect of damage persists.
f_plotFigure13(DSF_TAC, REF_FRAC)

Based on stiffness proxy

In a similar manner, the stiffness proxy is calculated for the earthquake signal. The same observations than for the TAC can be made: the values are lower than for low-amplitude vibrations - but the structure regains a higher stiffness when the amplitude of shaking increases. Unlike for the TAC, the values remain in the confidence interval - mostly due to the large scatter in the confidence interval that results from the low signal-to-noise ratio in the reference signals that undermine numerical intergration of acceleration to displacement.
Time = (0:1/fsamp:(length(WTr.DOF3)-1)/fsamp);
for Pairdummy = 1:length(TransmPairs)
for ModeDummy = 1:LMRef.nModes(Pairdummy)
sigtop = icwt(sgolayfilt(WTr.(['DOF',num2str(TransmPairs(Pairdummy,2))]), 3, 9), ...
fr_WTr, [FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1), FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2)]);
sigbot = icwt(sgolayfilt(WTr.(['DOF',num2str(TransmPairs(Pairdummy,1))]), 3, 9), ...
fr_WTr, [FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1), FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2)]);
DummyCounter = 1; IndexPairLow = find(Channels == TransmPairs(Pairdummy, 1));
for dummyLoop = IndexPairLow : nDOF
sigacc2(:, DummyCounter) = icwt(sgolayfilt(WTr.(['DOF',num2str(Channels(dummyLoop))]), 3, 9), ...
fr_WTr, [FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1), FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2)]);
DummyCounter = DummyCounter + 1;
end
sigaccFI = sum(sigacc2, 2); %Proxy of force (normalized by mass - which is considered unknown)
TopDisp = inte_frequency([Time', sigtop'], FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1), FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2), 2);
BaseDisp = inte_frequency([Time', sigbot'], FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1), FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2), 2);
RelDis = TopDisp - BaseDisp; %Proxy of relative displacement
[~, DSF_LFIT.(['Mode',num2str(ModeDummy)])(:,Pairdummy) , DSF_KPRX.(['Mode',num2str(ModeDummy)])(:,Pairdummy) , TimeVals, DSF_KPRX.Overall.(['Mode',num2str(ModeDummy)])(Pairdummy) ] = f_DSF_LC_batch(-1.*sigaccFI, RelDis', Time', fsamp, [1 35], WL, WL-OL);
end
end
f_plotFigure14(DSF_TAC, Kproxy, DSF_KPRX)

Spectral energy distribution

The same observation than for previous DSFs can also be made for the energy spread in the frequency domain.
for j = 1:length(TransmPairs)
for ModeDummy = 1:FSRef.nModes(j)
[~, M1EthOUT, ~] = f_WTindex_M1E_batch(WTr.(['DOF',num2str(TransmPairs(j,2))]), fr_WTr,LMRef.(['Pair',num2str(j)]).Peak(ModeDummy), WL, Time, (WL-OL));
[~, M1EthIN, tvals] = f_WTindex_M1E_batch(WTr.(['DOF',num2str(TransmPairs(j,1))]), fr_WTr,LMRef.(['Pair',num2str(j)]).Peak(ModeDummy), WL, Time, (WL-OL));
DSF_MER.Mean(j, ModeDummy) = mean(M1EthOUT ./ M1EthIN);
DSF_MER.StDev(j, ModeDummy) = std(M1EthOUT ./ M1EthIN);
DSF_MER.TimeHist.(['Mode',num2str(ModeDummy)])(:, j) = M1EthOUT ./ M1EthIN;
end
end
f_plotFigure16(DSF_TAC, DSF_MER, REF_MER)

Modal identification

The modal identification is done on the ambient signal (shown below) - to reduce the transient effects of the earthquake signal. Again, to facilitate picking the frequency peaks, the peaks are searched in a limited bandwidth of 0.6 to 1.2 times the reference frequency.
AMB_signal = WNM.Test3;
f_plotAmbientSig(AMB_signal)
WL = 30*fsamp; OL = 25*fsamp;
StartP = 1; EndP = WL; WindowDummy = 1;
while EndP <= length(AMB_signal)
[CPSD, fff] = f_CPSD_wind(AMB_signal(StartP:EndP, Channels), fsamp, 10*fsamp, 4*fsamp, 'hamming', 10*fsamp, 'c');
[Umat, Svec] = f_svd4fdd(CPSD, 1, 'act');
Svec = sgolayfilt(Svec, 3, 9);
EndIndex = f_findmatch(fff, 40);
% as we have predefined natural frequencies, we conduct a guided search in
% limited bands around these frequencies
for ModeDummy = 1: 2
ind1 = f_findmatch(fff, fnat_id(ModeDummy)*0.6);
ind2 = f_findmatch(fff, fnat_id(ModeDummy)*1.2);
[~, PeakIndexes(ModeDummy)] = max(Svec(ind1: ind2));
PeakIndexes(ModeDummy) = PeakIndexes(ModeDummy) + ind1 -1;
end
Ref_fnat.Ref2.fnat_id(WindowDummy, :) = fff(PeakIndexes);
MScompl = Umat(sort(PeakIndexes, 'ascend'), :, 1)';
% get the relative distance from the ground point (considered fix)
for j = 2:5
MSrel(j,1) = norm(MScompl(1,1)-MScompl(j,1)); MSrel(j,2) = norm(MScompl(1,2)-MScompl(j,2));
end
% get the relative angle from the ground point (considered fix)
MSangle = angle(MScompl); MSangle = MSangle - repmat(MSangle(1,:), 5, 1);
% construct and normalize the mode shape
Msha_sample = MSrel .* sign(cos(MSangle)); Msha_sample = Msha_sample ./ repmat(Msha_sample(5,:), 5, 1);
for ModeDummy = 1:2
for ChanDummy = 2:numel(unique(TransmPairs))-1
MSCinst.WN2.(['Mode',num2str(ModeDummy)])(WindowDummy, ChanDummy-1) = (Msha_sample(ChanDummy-1, ModeDummy) - 2*Msha_sample(ChanDummy, ModeDummy) + Msha_sample(ChanDummy+1, ModeDummy)) / (1.475*1.475);
end
MSCinst.WN2.(['Mode',num2str(ModeDummy)])(WindowDummy, numel(unique(TransmPairs))-1) = (Msha_sample(ChanDummy-1, ModeDummy) - Msha_sample(ChanDummy, ModeDummy) ) / (1.475*1.475);
end
StartP = StartP + (WL - OL); EndP = EndP + (WL - OL); WindowDummy = WindowDummy + 1;
end
% remove outliers but issue warning if too many outliers are detected
outliers = isoutlier(Ref_fnat.Ref2.fnat_id(:,1));
if (numel(find(outliers))/length(Ref_fnat.Ref2.fnat_id(:,1)))>0.15
disp('Warning! More than 15% of the data samples are identified as outliers')
end
f_plotFigure17(Ref_fnat, outliers)
As can be seen in the distribution of the natural frequencies that are identified, no significant change is observed in the first natural frequency, while a minor decrease of the second natural frequency can be observed. This observation can be quantidied with the Kullback-Leibler divergence that measured the distance between two datasets, which remain very low.
disp(['Kullback-Leibler divergence between healthy reference and new dataset: (Mode 1) ',num2str(f_KLdivergence(normpdf(linspace(4, 9, 101),mean(Ref_fnat.Ref0.fnat_id(:,1)), std(Ref_fnat.Ref0.fnat_id(:,1))), normpdf(linspace(4, 9, 101),mean(Ref_fnat.Ref2.fnat_id(:,1)), std(Ref_fnat.Ref2.fnat_id(:,1))))./100)])
Kullback-Leibler divergence between healthy reference and new dataset: (Mode 1) 0.12525
disp(['Kullback-Leibler divergence between healthy reference and new dataset: (Mode 2) ',num2str(f_KLdivergence(normpdf(linspace(19, 29, 101),mean(Ref_fnat.Ref0.fnat_id(:,2)), std(Ref_fnat.Ref0.fnat_id(:,2))), normpdf(linspace(19, 29, 101),mean(Ref_fnat.Ref2.fnat_id(:,2)), std(Ref_fnat.Ref2.fnat_id(:,2))))./100)])
Kullback-Leibler divergence between healthy reference and new dataset: (Mode 2) 0.6201
clearvars SigData WTr WTrAmp DSF_TAC DSF_CE DSF_KPRX DSF_LFIT DSF_MER DSF_MSC DSF_TAC_OverallMaha DSF_TS RefWeight

Near-real time damage assessment
The animation below shows the computation of DSFs during the earthquake.
The first row contains the ground motion acceleration and indicates whether an active earthquake is detected.
Rows 2-5 contain DSFs: TAC, KPRX, and MER - all evaluated for the sensor pair 0-4. Thresholds for the DSFs are defined based on the reference distribution. As discussed before, some nonlinearity and minor damage is observed during the earthquake - this is picked up by the DSFs falling out of the 95% confidence interval. However, at the end of the earthquake no residual damage is observed. In addition, the deviation during the earthquake remains small and does not affect all DSFs, therefore the green tag is attributed at the end.
The last row provides the multi-dimensional (square-root of the) Mahalonobis distance. It is evaluated for 6 DSFs: , , , , , and .In addition, based on the Mahalonobis distance, a tag is issued. Data-driven damage detection, in absence of data from damage classes, allows to split between healthy (green) and not healthy (amber or red), Hovever, based on the discrepancy from the reference dataset, amber or red tags may be issued as well.
ezgif.com-gif-maker.gif

3.2 Moderately strong earthquake (Test 5)

SigData = DATASET.Test5;
TimeVect = (1:length(SigData))./fsamp;
for ChanDummy = 1:nDOF
SigData(:, Channels(ChanDummy)) = f_butterworth(SigData(:, Channels(ChanDummy)), fsamp, 4, 35, 'low');
end
STA = 0.25; LTA = 5;
T1 = movavg(TimeVect', 'simple', floor(LTA*fsamp), 'fill');
T2 = movavg(TimeVect', 'simple', floor(STA*fsamp), 'fill');
MA1 = rms(SigData(:,1)) + movavg(SigData(:,1), 'simple', floor(LTA*fsamp), 'fill');
MA2 = rms(SigData(:,1)) + movavg(SigData(:,1), 'simple', floor(STA*fsamp), 'fill');
Ratio = abs(MA2) ./ abs(MA1);
Ratio = movmean(Ratio, 5);
Ind1 = min(find(Ratio < 0.7, 1, 'first'), find(Ratio > 1.5, 1, 'first')) - fsamp;
Ind2 = max(find(Ratio < 0.7, 1, 'last'), find(Ratio > 1.5, 1, 'last')) + 5*fsamp;
Ind2 = max(Ind2, Ind1+8*fsamp); % Impose a minimum length of 8 seconds
AnaRange = (Ind1 : Ind2);
f_plotFigure6(SigData,Channels,TimeVect,Ind1,Ind2,fsamp)
Sig2Analyse = SigData(Ind1: Ind2, :);
fb = cwtfilterbank('SignalLength',length(SigData),'SamplingFrequency',fsamp,...
'FrequencyLimits', fanaRange, 'VoicesPerOctave',48, 'TimeBandWidth',90);
for j = 1: nDOF
[wtrInstance, fr_WTr, CoI] = cwt(SigData(:, Channels(j)), 'FilterBank', fb);
WTr.(['DOF',num2str(Channels(j))]) = wtrInstance(:, Ind1:Ind2);
WTrAmp.(['DOF',num2str(Channels(j))]) = abs(wtrInstance(:, Ind1:Ind2));
end
fr_WTr_flud = flipud(fr_WTr);

Transmissibility

WL = 5*fsamp; OL = 4.5*fsamp;
StartP = 1; EndP = WL; WindowDummy = 1;
while EndP < length(Sig2Analyse)
for Pairdummy = 1:length(TransmPairs)
trans = sgolayfilt(WTrAmp.(['DOF',num2str(TransmPairs(Pairdummy,2))])(:, StartP:EndP), 3, 9) ...
./ sgolayfilt(WTrAmp.(['DOF',num2str(TransmPairs(Pairdummy,1))])(:, StartP:EndP), 3, 9);
RefRatio = 1 ./ sgolayfilt(WTrAmp.(['DOF',num2str(TransmPairs(Pairdummy,1))])(:, StartP:EndP), 3, 9); % this is used to avoid strong influence of the input spectrum
trans = flipud(trans);
[MScoh, fmsc] = mscohere(Sig2Analyse(StartP:EndP, TransmPairs(Pairdummy, 1)), Sig2Analyse(StartP:EndP, TransmPairs(Pairdummy, 2)), 2*fsamp, 1*fsamp, 10*fsamp, fsamp);
for ModeDummy = 1:TrRef.nModes(Pairdummy)
Index1 = f_findmatch(fr_WTr_flud, TrRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1));
Index2 = f_findmatch(fr_WTr_flud, TrRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2));
DSF_TAC.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy) = 1- f_MAC(...
TransmRef.(['Pair',num2str(Pairdummy)])(Index1:Index2), ...
mean(trans(Index1:Index2, :), 2));
% here a factor is introduced to indicate how much
% influence the input signal has on the transmissibility
FactorInput.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy) = f_MAC(...
mean(trans(Index1:Index2, :), 2), ...
mean(RefRatio(Index1:Index2,:),2));
DSF_TAC.mahadi.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy) = mahal(DSF_TAC.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy), REF_FRAC.Mode1(:, Pairdummy));
Index1 = f_findmatch(fmsc, TrRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1)); Index2 = f_findmatch(fmsc, TrRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2));
DSF_MSC.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy) = 1 - mean(MScoh(Index1:Index2));
DSF_MSC.mahadi.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy) = mahal(DSF_MSC.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy), REF_MSC.(['Mode',num2str(ModeDummy)])(:,Pairdummy));
end
end
StartP = StartP + (WL-OL); EndP = EndP + (WL-OL); WindowDummy = WindowDummy + 1;
end
Mahal.TAC.EQK5 = DSF_TAC.mahadi;
For localization purposes (see Section 4), an overall TAC is calculated and averaged over the entire length of the earthquake. To reduce the influence of the -highly variable and transient - input on the calculated transmissibility, a quality factor is calculated and used as a weighting factor when averaging the TAC over the entire seismic event:
,
where the MAC is the modal assurance criterion, denotes the amplitude of the spectral density of the frequency range of the signal and i refers to the input signal, while o designates the output signal. When the transmissibility has large peaks that result from valleys in the input frequency spectrum, the collinearity (i.e. the MAC value) will be high and the quality factor is small.
% we also compute here an overall value for the TAC - that will be used
% later on for localization
for Pairdummy = 1:length(TransmPairs)
DSF_TAC.OverallWeighted.(['Mode',num2str(ModeDummy)])(Pairdummy) = ...
sum(DSF_TAC.(['Mode',num2str(ModeDummy)])(:,Pairdummy) .* ...
(1-FactorInput.(['Mode',num2str(ModeDummy)])(:,Pairdummy))) / ...
sum((1-FactorInput.(['Mode',num2str(ModeDummy)])(:,Pairdummy)));
end
f_plotFigure13(DSF_TAC, REF_FRAC)
The second earthquake produces higher nonlinearity, when compared to the previous one. The TAC value excdeed twice the value that delimitates the 95%-confidence interval from healthy reference data. Therefore, onset of damage is concluded with high probability. In addition. although the TAC values decrease towars the end of the earthquake, when the amplitude of shaking decreases, it remains significantly above the threshold established with healthy reference data - indicating damage that permanently alters the structural condition.

Based on stiffness proxy

clearvars sigaccFI sigacc2 DSF_LFIT DSF_KPRX
Time = (0:1/fsamp:(length(WTr.DOF3)-1)/fsamp);
for Pairdummy = 1:length(TransmPairs)
for ModeDummy = 1:LMRef.nModes(Pairdummy)
sigtop = icwt(sgolayfilt(WTr.(['DOF',num2str(TransmPairs(Pairdummy,2))]), 3, 9), ...
fr_WTr, [FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1), FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2)]);
sigbot = icwt(sgolayfilt(WTr.(['DOF',num2str(TransmPairs(Pairdummy,1))]), 3, 9), ...
fr_WTr, [FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1), FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2)]);
DummyCounter = 1; IndexPairLow = find(Channels == TransmPairs(Pairdummy, 1));
for dummyLoop = IndexPairLow : nDOF
sigacc2(:, DummyCounter) = icwt(sgolayfilt(WTr.(['DOF',num2str(Channels(dummyLoop))]), 3, 9), ...
fr_WTr, [FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1), FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2)]);
DummyCounter = DummyCounter + 1;
end
sigaccFI = sum(sigacc2, 2); %Proxy of force (normalized by mass - which is considered unknown)
TopDisp = inte_frequency([Time', sigtop'], FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1), FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2), 2);
BaseDisp = inte_frequency([Time', sigbot'], FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1), FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2), 2);
RelDis = TopDisp - BaseDisp; %Proxy of relative displacement
[~, DSF_LFIT.(['Mode',num2str(ModeDummy)])(:,Pairdummy) , DSF_KPRX.(['Mode',num2str(ModeDummy)])(:,Pairdummy) , TimeVals, DSF_KPRX.Overall.(['Mode',num2str(ModeDummy)])(Pairdummy) ] = f_DSF_LC_batch(-1.*sigaccFI, RelDis', Time', fsamp, [1 35], WL, WL-OL);
end
end
f_plotFigure14(DSF_TAC, Kproxy, DSF_KPRX)
The KPRX feature confirms the the observation from the TAC: the structure undergoes significant nonlinearity and sustains damage.

Spectral energy distribution

for j = 1:length(TransmPairs)
for ModeDummy = 1:FSRef.nModes(j)
[~, M1EthOUT, ~] = f_WTindex_M1E_batch(WTr.(['DOF',num2str(TransmPairs(j,2))]), fr_WTr,LMRef.(['Pair',num2str(j)]).Peak(ModeDummy), WL, Time, (WL-OL));
[~, M1EthIN, tvals] = f_WTindex_M1E_batch(WTr.(['DOF',num2str(TransmPairs(j,1))]), fr_WTr,LMRef.(['Pair',num2str(j)]).Peak(ModeDummy), WL, Time, (WL-OL));
DSF_MER.Mean(j, ModeDummy) = mean(M1EthOUT ./ M1EthIN);
DSF_MER.StDev(j, ModeDummy) = std(M1EthOUT ./ M1EthIN);
DSF_MER.TimeHist.(['Mode',num2str(ModeDummy)])(:, j) = M1EthOUT ./ M1EthIN;
end
end
f_plotFigure16(DSF_TAC, DSF_MER, REF_MER)
Also the MER feature confirms significant nonlinearity and permanent damage.

Modal identification

Again, the ambient signal is used to compute the DSFs based on output-only modal identification.
AMB_signal = WNM.Test6;
f_plotAmbientSig(AMB_signal); xlim([0 90])
WL = 25*fsamp; OL = 20*fsamp;
StartP = 1; EndP = WL; WindowDummy = 1;
while EndP <= length(AMB_signal)
[CPSD, fff] = f_CPSD_wind(AMB_signal(StartP:EndP, Channels), fsamp, 10*fsamp, 7*fsamp, 'hamming', 10*fsamp, 'c');
[Umat, Svec] = f_svd4fdd(CPSD, 1, 'act');
Svec = sgolayfilt(Svec, 3, 9);
EndIndex = f_findmatch(fff, 40);
% as we have predefined natural frequencies, we conduct a guided search in
% limited bands around these frequencies
for ModeDummy = 1: 2
ind1 = f_findmatch(fff, fnat_id(ModeDummy)*0.4);
ind2 = f_findmatch(fff, fnat_id(ModeDummy)*1.1);
[~, PeakIndexes(ModeDummy)] = max(Svec(ind1: ind2));
PeakIndexes(ModeDummy) = PeakIndexes(ModeDummy) + ind1 -1;
end
Ref_fnat.Ref2.fnat_id(WindowDummy, :) = fff(PeakIndexes);
MScompl = Umat(sort(PeakIndexes, 'ascend'), :, 1)';
% get the relative distance from the ground point
for j = 2:5
MSrel(j,1) = norm(MScompl(1,1)-MScompl(j,1)); MSrel(j,2) = norm(MScompl(1,2)-MScompl(j,2));
end
% get the relative angle from the ground point
MSangle = angle(MScompl);
MSangle = MSangle - repmat(MSangle(1,:), 5, 1);
% construct the mode shape
Msha_sample = MSrel .* sign(cos(MSangle));
Msha_sample = Msha_sample ./ repmat(Msha_sample(5,:), 5, 1);
for ModeDummy = 1:2
for ChanDummy = 2:numel(unique(TransmPairs))-1
MSCinst.WN2.(['Mode',num2str(ModeDummy)])(WindowDummy, ChanDummy-1) = (Msha_sample(ChanDummy-1, ModeDummy) - 2*Msha_sample(ChanDummy, ModeDummy) + Msha_sample(ChanDummy+1, ModeDummy)) / (1.475*1.475);
end
MSCinst.WN2.(['Mode',num2str(ModeDummy)])(WindowDummy, numel(unique(TransmPairs))-1) = (Msha_sample(ChanDummy-1, ModeDummy) - Msha_sample(ChanDummy, ModeDummy) ) / (1.475*1.475);
end
StartP = StartP + (WL - OL); EndP = EndP + (WL - OL); WindowDummy = WindowDummy + 1;
end
outliers = isoutlier(Ref_fnat.Ref2.fnat_id(:,1));
if (numel(find(outliers))/length(Ref_fnat.Ref2.fnat_id(:,1)))>0.15
disp('Warning! More than 15% of the data samples are identified as outliers')
end
f_plotFigure17(Ref_fnat, outliers)
As can be seen in the figure, the fundamental frequencies are less sensitive to damage than the DSFs that are derived from the structure during the seismic event. However, two clusters form for the first natural frequency, indicating a deviation from the healthy reference state. In addition, the second frequency deviates from the reference distribution, as also confirmed by the higher Kullback-Leibler divergence as computed below.
f_KLdivergence(normpdf(linspace(4, 9, 101),mean(Ref_fnat.Ref0.fnat_id(:,1)), std(Ref_fnat.Ref0.fnat_id(:,1))), normpdf(linspace(4, 9, 101),mean(Ref_fnat.Ref2.fnat_id(:,1)), std(Ref_fnat.Ref2.fnat_id(:,1))))./100
ans = 0.2943
f_KLdivergence(normpdf(linspace(20, 30, 101),mean(Ref_fnat.Ref0.fnat_id(:,2)), std(Ref_fnat.Ref0.fnat_id(:,2))), normpdf(linspace(20, 30, 101),mean(Ref_fnat.Ref2.fnat_id(:,2)), std(Ref_fnat.Ref2.fnat_id(:,2))))./100
ans = 0.4528
Near-real time damage detection
The derivation of DSFs over small time windows allows to track the structural state in near-real time. For this earthquake example, the data indicates damage (which is in agreement with the observation on the test specimen). Again, it is noted that for data-driven anomaly detection, the amber/red tags may give some indication about the severity of damage, but are mostly correlated with the certainty, with which the damage is detected. In this case, the red tag (of the bottom subplot) indicates that - based on the DSFs - damage is detected with high confidence.
Earthquake5_Demo.gif
% store some data for the localization section and clear the others.
DSF2.KPRX = DSF_KPRX; DSF2.TAC = DSF_TAC.OverallWeighted; DSF2.MSC = MSCinst.WN2;
clearvars SigData WTr WTrAmp DSF_TAC DSF_CE DSF_KPRX DSF_LFIT DSF_MER DSF_MSC DSF_TAC_OverallMaha DSF_TS RefWeight

3.3 Strong earthquake

SigData = DATASET.Test8;
TimeVect = (1:length(SigData))./fsamp;
for ChanDummy = 1:nDOF
SigData(:, Channels(ChanDummy)) = f_butterworth(SigData(:, Channels(ChanDummy)), fsamp, 4, 35, 'low');
end
STA = 0.25; LTA = 5;
T1 = movavg(TimeVect', 'simple', floor(LTA*fsamp), 'fill');
T2 = movavg(TimeVect', 'simple', floor(STA*fsamp), 'fill');
MA1 = rms(SigData(:,1)) + movavg(SigData(:,1), 'simple', floor(LTA*fsamp), 'fill');
MA2 = rms(SigData(:,1)) + movavg(SigData(:,1), 'simple', floor(STA*fsamp), 'fill');
Ratio = abs(MA2) ./ abs(MA1);
Ratio = movmean(Ratio, 5);
Ind1 = min(find(Ratio < 0.7, 1, 'first'), find(Ratio > 1.5, 1, 'first')) - fsamp;
Ind2 = max(find(Ratio < 0.7, 1, 'last'), find(Ratio > 1.5, 1, 'last')) + 5*fsamp;
Ind2 = max(Ind2, Ind1+8*fsamp); % Impose a minimum length of 8 seconds
AnaRange = (Ind1 : Ind2);
f_plotFigure6(SigData,Channels,TimeVect,Ind1,Ind2,fsamp)
Sig2Analyse = SigData(Ind1: Ind2, :);
fb = cwtfilterbank('SignalLength',length(SigData),'SamplingFrequency',fsamp,...
'FrequencyLimits', fanaRange, 'VoicesPerOctave',48, 'TimeBandWidth',90);
for j = 1: nDOF
[wtrInstance, fr_WTr, CoI] = cwt(SigData(:, Channels(j)), 'FilterBank', fb);
WTr.(['DOF',num2str(Channels(j))]) = wtrInstance(:, Ind1:Ind2);
WTrAmp.(['DOF',num2str(Channels(j))]) = abs(wtrInstance(:, Ind1:Ind2));
end
fr_WTr_flud = flipud(fr_WTr);

Transmissibility

WL = 5*fsamp; OL = 4.5*fsamp;
StartP = 1; EndP = WL; WindowDummy = 1;
while EndP < length(Sig2Analyse)
for Pairdummy = 1:length(TransmPairs)
trans = sgolayfilt(WTrAmp.(['DOF',num2str(TransmPairs(Pairdummy,2))])(:, StartP:EndP), 3, 9) ...
./ sgolayfilt(WTrAmp.(['DOF',num2str(TransmPairs(Pairdummy,1))])(:, StartP:EndP), 3, 9);
RefRatio = 1 ./ sgolayfilt(WTrAmp.(['DOF',num2str(TransmPairs(Pairdummy,1))])(:, StartP:EndP), 3, 9); % this is used to avoid strong influence of the input spectrum
trans = flipud(trans);
[MScoh, fmsc] = mscohere(Sig2Analyse(StartP:EndP, TransmPairs(Pairdummy, 1)), Sig2Analyse(StartP:EndP, TransmPairs(Pairdummy, 2)), 2*fsamp, 1*fsamp, 10*fsamp, fsamp);
for ModeDummy = 1:TrRef.nModes(Pairdummy)
Index1 = f_findmatch(fr_WTr_flud, TrRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1));
Index2 = f_findmatch(fr_WTr_flud, TrRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2));
DSF_TAC.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy) = 1- f_MAC(...
TransmRef.(['Pair',num2str(Pairdummy)])(Index1:Index2), ...
mean(trans(Index1:Index2, :), 2));
DSF_TAC.mahadi.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy) = mahal(DSF_TAC.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy), REF_FRAC.Mode1(:, Pairdummy));
Index1 = f_findmatch(fmsc, TrRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1)); Index2 = f_findmatch(fmsc, TrRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2));
DSF_MSC.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy) = 1 - mean(MScoh(Index1:Index2));
DSF_MSC.mahadi.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy) = mahal(DSF_MSC.(['Mode',num2str(ModeDummy)])(WindowDummy,Pairdummy), REF_MSC.(['Mode',num2str(ModeDummy)])(:,Pairdummy));
end
end
StartP = StartP + (WL-OL); EndP = EndP + (WL-OL); WindowDummy = WindowDummy + 1;
end
Mahal.TAC.EQK8 = DSF_TAC.mahadi;
f_plotFigure13(DSF_TAC, REF_FRAC)
The distance between the 95%-confidence interval from healthy reference data and the DSF computed from the earthquake signal is further increased, when compared with the previous case. This indicates more severe damage that is very likely to have an influence on the residual capacity of the building.

Based on stiffness proxy

clearvars sigaccFI sigacc2 DSF_LFIT DSF_KPRX
Time = (0:1/fsamp:(length(WTr.DOF3)-1)/fsamp);
for Pairdummy = 1:length(TransmPairs)
for ModeDummy = 1:LMRef.nModes(Pairdummy)
sigtop = icwt(sgolayfilt(WTr.(['DOF',num2str(TransmPairs(Pairdummy,2))]), 3, 9), ...
fr_WTr, [FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1), FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2)]);
sigbot = icwt(sgolayfilt(WTr.(['DOF',num2str(TransmPairs(Pairdummy,1))]), 3, 9), ...
fr_WTr, [FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1), FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2)]);
DummyCounter = 1; IndexPairLow = find(Channels == TransmPairs(Pairdummy, 1));
for dummyLoop = IndexPairLow : nDOF
sigacc2(:, DummyCounter) = icwt(sgolayfilt(WTr.(['DOF',num2str(Channels(dummyLoop))]), 3, 9), ...
fr_WTr, [FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1), FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2)]);
DummyCounter = DummyCounter + 1;
end
sigaccFI = sum(sigacc2, 2); %Proxy of force (normalized by mass - which is considered unknown)
TopDisp = inte_frequency([Time', sigtop'], FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1), FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2), 2);
BaseDisp = inte_frequency([Time', sigbot'], FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 1), FSRef.(['Pair',num2str(Pairdummy)]).fRange(ModeDummy, 2), 2);
RelDis = TopDisp - BaseDisp; %Proxy of relative displacement
[~, DSF_LFIT.(['Mode',num2str(ModeDummy)])(:,Pairdummy) , DSF_KPRX.(['Mode',num2str(ModeDummy)])(:,Pairdummy) , TimeVals, DSF_KPRX.Overall.(['Mode',num2str(ModeDummy)])(Pairdummy) ] = f_DSF_LC_batch(-1.*sigaccFI, RelDis', Time', fsamp, [1 35], WL, WL-OL);
end
end
f_plotFigure14(DSF_TAC, Kproxy, DSF_KPRX)
The stiffness approximator (KPRX) confirms the findings from the TAC and shows stiffness values in the range of 0 - indicating very strong nonlinearity and important plastic deformations. Again, this underlines the likely reduction of seismic capacity.

Spectral energy distribution

for j = 1:length(TransmPairs)
for ModeDummy = 1:FSRef.nModes(j)
[~, M1EthOUT, ~] = f_WTindex_M1E_batch(WTr.(['DOF',num2str(TransmPairs(j,2))]), fr_WTr,LMRef.(['Pair',num2str(j)]).Peak(ModeDummy), WL, Time, (WL-OL));
[~, M1EthIN, tvals] = f_WTindex_M1E_batch(WTr.(['DOF',num2str(TransmPairs(j,1))]), fr_WTr,LMRef.(['Pair',num2str(j)]).Peak(ModeDummy), WL, Time, (WL-OL));
DSF_MER.Mean(j, ModeDummy) = mean(M1EthOUT ./ M1EthIN);
DSF_MER.StDev(j, ModeDummy) = std(M1EthOUT ./ M1EthIN);
DSF_MER.TimeHist.(['Mode',num2str(ModeDummy)])(:, j) = M1EthOUT ./ M1EthIN;
end
end
f_plotFigure16(DSF_TAC, DSF_MER, REF_MER)
Also the MER shows stronger nonlinearity and thus, damage, than for the precious earthquake.

Modal identification

AMB_signal = WNM.Test9;
WL = 25*fsamp; OL = 20*fsamp;
StartP = 1; EndP = WL; WindowDummy = 1;
while EndP <= length(AMB_signal)
[CPSD, fff] = f_CPSD_wind(AMB_signal(StartP:EndP, Channels), fsamp, 10*fsamp, 7.5*fsamp, 'hamming', 10*fsamp, 'c');
[Umat, Svec] = f_svd4fdd(CPSD, 1, 'act');
Svec = sgolayfilt(Svec, 3, 9);
EndIndex = f_findmatch(fff, 40);
% as we have predefined natural frequencies, we conduct a guided search in
% limited bands around these frequencies
for ModeDummy = 1: 2
ind1 = f_findmatch(fff, fnat_id(ModeDummy)*0.30);
ind2 = f_findmatch(fff, fnat_id(ModeDummy)*1.05);
[~, PeakIndexes(ModeDummy)] = max(Svec(ind1: ind2));
PeakIndexes(ModeDummy) = PeakIndexes(ModeDummy) + ind1 -1;
end
Ref_fnat.Ref2.fnat_id(WindowDummy, :) = fff(PeakIndexes);
MScompl = Umat(sort(PeakIndexes, 'ascend'), :, 1)';
% get the relative distance from the ground point
for j = 2:5
MSrel(j,1) = norm(MScompl(1,1)-MScompl(j,1)); MSrel(j,2) = norm(MScompl(1,2)-MScompl(j,2));
end
% get the relative angle from the ground point
MSangle = angle(MScompl);
MSangle = MSangle - repmat(MSangle(1,:), 5, 1);
% construct the mode shape
Msha_sample = MSrel .* sign(cos(MSangle));
Msha_sample = Msha_sample ./ repmat(Msha_sample(5,:), 5, 1);
for ModeDummy = 1:2
for ChanDummy = 2:numel(unique(TransmPairs))-1
MSCinst.WN2.(['Mode',num2str(ModeDummy)])(WindowDummy, ChanDummy-1) = (Msha_sample(ChanDummy-1, ModeDummy) - 2*Msha_sample(ChanDummy, ModeDummy) + Msha_sample(ChanDummy+1, ModeDummy)) / (1.475*1.475);
end
MSCinst.WN2.(['Mode',num2str(ModeDummy)])(WindowDummy, numel(unique(TransmPairs))-1) = (Msha_sample(ChanDummy-1, ModeDummy) - Msha_sample(ChanDummy, ModeDummy) ) / (1.475*1.475);
end
StartP = StartP + (WL - OL); EndP = EndP + (WL - OL); WindowDummy = WindowDummy + 1;
end
outliers = isoutlier(Ref_fnat.Ref2.fnat_id(:,1));
if (numel(find(outliers))/length(Ref_fnat.Ref2.fnat_id(:,1)))>0.15
disp('Warning! More than 15% of the data samples are identified as outliers')
end
f_plotFigure17(Ref_fnat, outliers)
In this case, the damage is very severe and significantly alters the natural frequencies of the building, leading to a frequency reduction of almost 50% for the first natural frequency. The strong discrepancy between the healthy and the damaged distribution can also be shown through the Kullback-Leibler divergence, computed below.
f_KLdivergence(normpdf(linspace(4, 9, 101),mean(Ref_fnat.Ref0.fnat_id(:,1)), std(Ref_fnat.Ref0.fnat_id(:,1))), normpdf(linspace(4, 9, 101),mean(Ref_fnat.Ref2.fnat_id(:,1)), std(Ref_fnat.Ref2.fnat_id(:,1))))./100
ans = 7.5604
f_KLdivergence(normpdf(linspace(15, 30, 201),mean(Ref_fnat.Ref0.fnat_id(:,2)), std(Ref_fnat.Ref0.fnat_id(:,2))), normpdf(linspace(15, 30, 201),mean(Ref_fnat.Ref2.fnat_id(:,2)), std(Ref_fnat.Ref2.fnat_id(:,2))))./200
ans = 6.5465
clearvars SigData WTr WTrAmp DSF_TAC DSF_CE DSF_KPRX DSF_LFIT DSF_MER DSF_MSC DSF_TAC_OverallMaha DSF_TS RefWeight

3.4 Summary of data-driven damage detection

When merging the information from the three analysed earthquakes, the progression of damage can be seen; as shown in the figure below that tracks the Mahalonobis distance between the and the corresponding healthy reference distribution. The structure undergoes stronger nonlinearity during the largest amplitudes of he earthquake and thus, the feature deviates stronger from the reference state. When the amplitude of shaking reduce - towards the end of each earthquake - the distance between the damaged state and the healthy state also slightly reduces. This is most likely due to the amplitude-dependency of masonry (and reinforced concrete) structures and the stiffness recovery due to crack closure after the earthquake.
f_plotFigure19(Mahal)

4. Data-driven Damage Localization

For the case of a moderate earthquake (see Section 2.3), the capabilities of damage-sensitive features to locate the damage is tested.
First, the stiffness proxy can be used to compute the change with respect to the reference value. Although the stiffness proxy is very noisy - and the relative displacement on the last floor is too small to compute relevant values - this DSF indicates that most damage is on the first floor and that from the third floor upwards no damage is expected. This again corresponds to the reality that has been observed on the specimen.
f_plotLocal1(Kproxy, DSF2)
Similarely, the TAC indicates that the most severe damage is to be expected in the first 2 floors - which is in agreement with the observations on the test specimen.
f_plotLocal2(DSF2)

5. Damage quantification

5.1 Physics-based nonlinear simulation model

Damage quantification requires a physics-based model to pre-establish a link between the value of damage-sensitive features and the probability of exceeding discrete damage states, which are commonly used in regional loss assessment to quantify earthquake damage. In this application, a simplified two-dimensional MDOF model is used to construct sucht fragility functions based on DSFs. A schematical representation of the model is shown below.
Model.png
An example of a fragility curve based on a damage-sensitive feature (DSF) - in this case the TAC - is shown below. Three damage states are defined and correspond to warning states:
The near-real time tagging of buildings as either undamaged or safe for occupancy (green) may reduce the downtime of buildings that have not been inspected yet and therefore, increase the resilience of communities with respect to earthquakes.
Using the DSFs that have been computed in Section 3, the following tags are predicted:
The model-based damage quantification for earthquake 2 confirms the absence of significant (residual) damage that has been concluded with the data-driven approach in section 3.1. In addition, the modal allows to refine the tagging of buildings, for which the presence of damage has been detected, into moderate and severe damage. This helps in post-earthquake decision-making and is a starting step to quantify the residual seismic performance of the damaged building.

References

Beyer, K., Tondelli, M., Petry, S. and Peloso, S., 2015. "Dynamic testing of a four-storey building with reinforced concrete and unreinforced masonry walls: prediction, test results and data set." Bulletin of Earthquake Engineering, 13(10), pp.3015-3064.
Mahalonobis, P.C.. 1936. "On the generalised distance in statistics" Proceedings of the National Institute of Science of India, 2, 49-55.
Martakis, P., Reuland, Y. and Chatzi, E., 2021. "Amplitude Dependency Effects in the Structural Identification of Historic Masonry Buildings". In International Conference of the European Association on Quality Control of Bridges and Structures (pp. 140-147). Springer.
Noh, H.Y., Nair, K.K., Lignos, D.G. and Kiremidjian, A.S., 2011. Use of wavelet-based damage-sensitive features for structural damage diagnosis using strong motion data. Journal of Structural Engineering, 137(10), pp.1215-1228.
Reuland, Y., Martakis, P., Chatzi, E., 2022. "Deliverable 4.5: The use of structural health monitoring for rapid loss assessment" Horizon2020 project RISE.
Rytter, A., 1993, “Vibration based inspection of civil engineering structures,” PhD Dissertation, Department of Building Technology and Structural Engineering, AalborgUniversity, Denmark.
Tondelli, Marco, Petry, Sarah, Peloso, Simone, & Beyer, Katrin. (2014). "Dynamic testing of a four-storey building with reinforced concrete and unreinforced masonry wall: Data set [Data set]." https://doi.org/10.5281/zenodo.11578
Zhou, Y.L., Figueiredo, E., Maia, N. and Perera, R., 2015. "Damage detection and quantification using transmissibility coherence analysis." Shock and Vibration, 2015.